fix(solo): MultiGeneUMI_CR decides ownership on corrected UMIs (stacked on #174) - #175
fix(solo): MultiGeneUMI_CR decides ownership on corrected UMIs (stacked on #174)#175BenjaminDEMAILLE wants to merge 4 commits into
Conversation
`--soloUMIfiltering MultiGeneUMI_CR` kept every gene tied at the highest
read count. CellRanger's rule is the opposite on exactly that case: the
gene with the *strictly* highest count takes the UMI, and a tie means no
gene counts it.
STAR walks the genes keeping a running maximum and clears its winner
whenever it meets an equal count
(`SoloFeature_collapseUMIall.cpp:212-224`):
if (ig.second>maxu) { maxu=ig.second; maxg=ig.first; }
else if (ig.second==maxu) { maxg=-1; };
...
if ( maxg+1==0 ) continue; // not counted for any gene
One read per gene is the ordinary shape of a multi-gene UMI, and it is
always a tie, so the old rule made the flag inert in practice rather
than merely inaccurate. Measured on a 20 000-read 10x fixture (200 cells
from the real v3 whitelist, 400 genes, 720 UMIs deliberately shared
between two genes), against STAR 2.7.11b with the same flags:
identical entries STAR counts rustar counts
before 13 749 / 14 806 15 423 16 465
after 13 902 / 13 967 15 423 15 414
The flag removed nothing at all before; STAR removes 1 030 counts. The
gap goes from +1 042 to -9.
The outcome does not depend on the order the genes are visited — a
strict maximum always ends as the winner, a tie always ends with none —
so iterating a `HashMap` here stays deterministic.
`multi_gene_umi_cr_drops_a_tie_entirely` pins the case the old tests
missed: they only covered 3 reads against 1, where both rules agree.
Not yet implemented, and stated so rather than left to be discovered:
STAR applies a second condition, that the winning gene must also hold
the top count among *uncorrected* UMIs (`umiGeneMapCount0`, same file,
lines 226-232). That needs the pre-correction counts, which this code
does not keep. The 65 entries still differing out of 13 967 are the
place to look for its effect.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…w matrix
STARsolo's raw matrix has a column per whitelist barcode. For 10x v3 that
is 3 686 400 columns and a 62 MB `barcodes.tsv`, nearly all zeros.
CellRanger's `raw_feature_bc_matrix` has a column per *observed* barcode.
The two files therefore share no keys, which is not a rounding difference
in a comparison, it is zero overlap: comparing our raw output against a
real `cellranger count` run gave 0 identical entries out of 27 396 until
the columns were reconciled.
`--soloOutRawBarcodes Observed` narrows the raw matrix to the barcodes
that carry a count. Default `Whitelist` keeps what STARsolo writes, so
nothing changes for anyone not asking.
Measured on the 20 000-read fixture:
Whitelist 3 686 400 barcodes barcodes.tsv 62 668 800 bytes
Observed 200 barcodes barcodes.tsv 3 400 bytes
with identical counts on both sides: 13 937 entries, 15 414 counts.
`finalize_matrix` already took a column remap for the filtered matrix, so
this reuses it rather than adding a second path. The observed set is read
back from the streamed body, which costs one pass and only when the flag
is on.
This is a **non-STAR flag** and needs sign-off; recorded in
`DIVERGENCE.md` §3.2 rather than presented as parity.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
STAR corrects UMIs within each gene *before* deciding which gene owns a UMI, and applies two conditions, not one (`SoloFeature_collapseUMIall.cpp:134-148` and `:203-235`): 1. one gene must hold a strictly higher read count than every other, on the **corrected** UMI map — that is scverse#173, already landed; 2. and that winner must not be beaten in the **uncorrected** map at the same key. The second condition exists because correction moves reads between UMIs: a gene can win only because correction folded a neighbouring UMI onto it, and STAR rejects that win rather than counting it. Reproducing it needs the order STAR uses. The generic path here filters multi-gene UMIs first and corrects afterwards, which cannot express either condition: by the time correction happens the ownership decision is already made. `MultiGeneUMI_CR` therefore takes its own path, which is also what STAR does — the flag is only valid with `--soloUMIdedup 1MM_CR`, so there is no combination this bypasses. `cellranger_1mm_map` exposes the correction mapping that `cellranger_1mm` already computed and threw away. Measured against **CellRanger 10.0.0** on the 20 000-read fixture from scverse#172, with scverse#165 and scverse#173 also applied: identical entries CellRanger rustar scverse#165 + scverse#173 13 651 / 13 709 15 111 15 091 plus this change 13 676 / 13 709 15 111 15 116 Entries CellRanger has and we do not go from 29 to 7, and the count gap from -20 to +5, which is 0.03%. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Correction: the CellRanger numbers in this PR do not reproduceI re-derived the whole comparison from a clean state today, both sides in one pass, and the absolute figures I published are wrong. The relative effects hold. The headline claim does not, and it was the headline, so this needs saying plainly rather than in a footnote. Everything below is from today, one fixture, one
What survives
What does not survive"within 0.03%", "15 116", and "rustar is closer to CellRanger than STAR is". All three are wrong. rustar is at +2.17%; STAR is at +0.09%. Compared directly against STAR under identical flags, rustar is +315 counts with 273 entries STAR does not have. That is a real open divergence, and I reported it as closed. How it happenedThe fixture reference was regenerated between that measurement and this one, and I compared a rustar run against a CellRanger run derived from a different state of it. The relative deltas were unaffected, which is why the per-PR effects still reproduce and I did not notice. I should have re-derived both sides in one pass before publishing an absolute figure, and from here I will. What this changes about the PRNothing about the code, and nothing about the argument for it: the flags still move the matrix from 8.96% to 2.17% off CellRanger, which is most of the gap and is what the change is for. What changes is that the remaining 2.17% is open, not closed, and the 273 entries rustar has that STAR does not are the next thing to chase. I would rather that be visible in the review than discovered later. I have edited the numbers in the PR description to match this table. |
Retracting the correction aboveThe correction I posted earlier was itself wrong, and I would rather say so immediately than leave it standing. The original numbers in this PR were right. What I failed to hold constant when I re-measured was #165, which is an open PR and not in this stack's base. Re-measured today, both sides in one pass, one fixture, one
The last row reproduces the original claim exactly, including the 13 676 / 13 709 identical entries. The PR body already named the precondition — it said "with #165 and #173" — and I dropped that condition when re-measuring, then published the resulting worse number as a correction of a claim that was never wrong. So, to be unambiguous about what holds:
Practical consequence for review: #165 should merge before this stack, or the numbers in these PRs will not reproduce. I have restored the original figures in the description and made that dependency explicit rather than parenthetical. The lesson I am taking from it, since it cost you reading time twice: when a measurement disagrees with a published one, the first thing to check is what moved between them, not the published number. |
* fix(solo): MultiGeneUMI_CR gives a tied UMI to nobody, not to everybody
`--soloUMIfiltering MultiGeneUMI_CR` kept every gene tied at the highest
read count. CellRanger's rule is the opposite on exactly that case: the
gene with the *strictly* highest count takes the UMI, and a tie means no
gene counts it.
STAR walks the genes keeping a running maximum and clears its winner
whenever it meets an equal count
(`SoloFeature_collapseUMIall.cpp:212-224`):
if (ig.second>maxu) { maxu=ig.second; maxg=ig.first; }
else if (ig.second==maxu) { maxg=-1; };
...
if ( maxg+1==0 ) continue; // not counted for any gene
One read per gene is the ordinary shape of a multi-gene UMI, and it is
always a tie, so the old rule made the flag inert in practice rather
than merely inaccurate. Measured on a 20 000-read 10x fixture (200 cells
from the real v3 whitelist, 400 genes, 720 UMIs deliberately shared
between two genes), against STAR 2.7.11b with the same flags:
identical entries STAR counts rustar counts
before 13 749 / 14 806 15 423 16 465
after 13 902 / 13 967 15 423 15 414
The flag removed nothing at all before; STAR removes 1 030 counts. The
gap goes from +1 042 to -9.
The outcome does not depend on the order the genes are visited — a
strict maximum always ends as the winner, a tie always ends with none —
so iterating a `HashMap` here stays deterministic.
`multi_gene_umi_cr_drops_a_tie_entirely` pins the case the old tests
missed: they only covered 3 reads against 1, where both rules agree.
Not yet implemented, and stated so rather than left to be discovered:
STAR applies a second condition, that the winning gene must also hold
the top count among *uncorrected* UMIs (`umiGeneMapCount0`, same file,
lines 226-232). That needs the pre-correction counts, which this code
does not keep. The 65 entries still differing out of 13 967 are the
place to look for its effect.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
* docs(changelog): record the MultiGeneUMI_CR tie fix
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
* fix(solo): MultiGeneUMI_CR decides ownership on corrected UMIs
STAR corrects UMIs within each gene *before* deciding which gene owns a
UMI, and applies two conditions, not one
(`SoloFeature_collapseUMIall.cpp:134-148` and `:203-235`):
1. one gene must hold a strictly higher read count than every other, on
the **corrected** UMI map — that is #173, already landed;
2. and that winner must not be beaten in the **uncorrected** map at the
same key.
The second condition exists because correction moves reads between UMIs:
a gene can win only because correction folded a neighbouring UMI onto it,
and STAR rejects that win rather than counting it.
Reproducing it needs the order STAR uses. The generic path here filters
multi-gene UMIs first and corrects afterwards, which cannot express either
condition: by the time correction happens the ownership decision is
already made. `MultiGeneUMI_CR` therefore takes its own path, which is
also what STAR does — the flag is only valid with `--soloUMIdedup 1MM_CR`,
so there is no combination this bypasses.
`cellranger_1mm_map` exposes the correction mapping that
`cellranger_1mm` already computed and threw away.
Measured against **CellRanger 10.0.0** on the 20 000-read fixture from
#172, with #165 and #173 also applied:
identical entries CellRanger rustar
#165 + #173 13 651 / 13 709 15 111 15 091
plus this change 13 676 / 13 709 15 111 15 116
Entries CellRanger has and we do not go from 29 to 7, and the count gap
from -20 to +5, which is 0.03%.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
* docs(divergence): drop a reference to a test that no longer exists
* feat(solo): bit-exact libc++ mt19937, generate_canonical and discrete_distribution
STARsolo's `EmptyDrops_CR` rescue draws from `std::mt19937`, converts to
doubles with `std::generate_canonical<double, 53>`, and picks categories with
`std::discrete_distribution`. Two of those three are implementation-defined in
the parts that matter: the standard fixes mt19937's output but not how
`generate_canonical` consumes it, and says nothing about how
`discrete_distribution` maps a uniform onto categories.
So porting "the algorithm" is not enough — it has to be libc++'s algorithm,
because that is what STAR is built against and where its numbers come from.
libc++ accumulates two 32-bit draws in *ascending* significance and divides by
2^64; a most-significant-first accumulation, or one draw scaled to 53 bits,
both give perfectly good uniforms and neither reproduces STAR.
Every expected value in the tests came out of a C++ program compiled against
the real libc++ and run, not from reading its source. `tests/libcxx_oracle.cpp`
is that program, kept so the values can be regenerated rather than trusted.
`generate_canonical` is compared as bit patterns, since a difference in the
last place changes which category a sample lands in.
Not yet wired into the EmptyDrops path. `solo::count` samples with a
`SplitMix64` stream under a comment calling it "WeightedIndex-equivalent;
empirically byte-identical EmptyDrops cell calls" — a claim that cannot hold in
general, since two unrelated generators cannot agree on an arbitrary number of
draws. It is true of whatever was checked and unknown elsewhere. Replacing it
moves cell calls, so it belongs in its own change with the solo differential
run against it.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
* feat(solo): EmptyDrops_CR uses Simple Good-Turing and libc++'s sampler
Two approximations in the CellRanger cell-calling path are replaced by what
CellRanger and STAR actually compute. Both move cell calls, which is the point:
the previous numbers were plausible rather than right.
The ambient profile is now smoothed with Simple Good-Turing (Gadsby & Sampson,
via Elworthy's implementation, which is what STAR vendors). The ambient counts
come from a small sample of empty droplets, so a gene seen twice there is not
twice as likely as one seen once, and a gene seen zero times is not impossible —
it is one the sample was too small to show. SGT fits the frequency spectrum and
reserves mass for the unseen from the singleton rate, then smooths the rest
along a log-log line. What was here before had the right shape and the wrong
numbers: it reserved mass the same way but distributed the remainder in
proportion to raw counts, with no smoothing at all.
The Monte-Carlo null is now drawn with libc++'s `std::mt19937` and
`std::discrete_distribution`, seeded `19760110 * (isim + 1)` per simulation, as
STAR seeds it. The previous sampler was a SplitMix64 stream under a comment
calling it "WeightedIndex-equivalent; empirically byte-identical EmptyDrops cell
calls" — a claim that cannot hold in general, since two unrelated generators
cannot agree over an arbitrary number of draws. The libc++ types were ported and
checked against real libc++ in the previous commit on this branch; this wires
them in. One generator per simulation, no shared state, so the walks still run
in any order on any number of threads and give the same p-values.
D17 comes with it: STAR leaves `PZero` uninitialised when the spectrum has fewer
than five distinct frequencies and `analyse()` bails, so it reads whatever the
stack held. Here it is zero from construction, which is what "no basis for
reserving unseen mass" means. Recorded in docs-old/dev/divergences.md.
* docs: record the EmptyDrops SGT divergence in DIVERGENCE.md
Section 1.2, in the What STAR does / What rustar-aligner does / Why / Impact /
Source format CONTRIBUTING.md asks for, replacing the docs-old file the earlier
version of this work carried.
* docs(divergence): file the EmptyDrops entry under section 1, note the
second RNG
* fix(params): refuse MultiGeneUMI_CR without --soloUMIdedup 1MM_CR
* feat(solo): --soloFeatures Transcript3p, with --soloClusterCBfile
Quantifies transcripts rather than genes, from where each read's 3' end sits
relative to each transcript's. In a 3'-biased assay that distance is what
separates isoforms: a read 200 bases from the end of one and 4000 from the end
of another is evidence for the first. The distribution of those distances is
estimated from the run's own histogram, smoothed and cut where the 3' peak
decays into the body, and used as the likelihood in an EM over UMIs.
Concordance needed no new code. `align_to_transcripts` already refuses to
project an alignment that leaves the transcript, touches an intron, or crosses a
junction the transcript does not have — which is exactly STAR's `Concordant`
(`Transcriptome_classifyAlign.cpp`). A projection that survives is concordant;
one that does not, is not. The projection also puts the 5' end at coordinate
zero for both strands, so the distance to the 3' end is one expression rather
than two.
Two behaviours worth stating because they are not the obvious ones:
Output is per cluster, not per cell, and `--soloClusterCBfile` is required.
A single cell does not have enough UMIs to resolve isoforms, so the EM would be
fitting noise. Asking for the feature without a clustering is refused rather
than run.
A UMI seen on several reads contributes the *intersection* of their transcript
sets. Those reads came from one molecule, so a transcript missing from any of
them cannot be its source. Taking the union would let a single stray read
resurrect an isoform every other read excluded.
Two of STAR's quirks are reproduced rather than corrected, because the cut point
and every weight depend on them: the running-average divisor is `min(2N+1,
i + N)` rather than the number of elements actually summed, and the transcript
length factor is taken from the cumulative distribution at `trLen - 1`
(`SoloFeature_quantTranscript.cpp`).
Numbers are formatted the way C++'s default stream prints them — six
significant digits, fixed inside `[1e-4, 1e6)` and scientific outside — since
the normalised distribution runs down to ~1e-4 where Rust's `{}` and C++'s
default disagree on both notation and digit count.
* refactor(solo): drop Transcript3pAcc::merge, which nothing calls
Records are accumulated under a mutex, so there are no partials to merge. It was
dead from the moment it was written; CONTRIBUTING.md rules out shipping it.
* docs(solo): note that STAR marks Transcript3p under development
parametersDefault puts both Transcript3p and --soloClusterCBfile between
"#####UnderDevelopment_begin : not supported - do not use" and
"#####UnderDevelopment_end", and STAR --help prints that banner around
them. The module said none of this.
It matters for how the port is read: it follows STAR's code, so it
inherits the unfinished parts of that code, and a differential against
STAR compares two implementations of something STAR does not support.
A reviewer should be told that before deciding to take it.
* fix(solo): implement MultiGeneUMI_All instead of aliasing it to MultiGeneUMI
`--soloUMIfiltering MultiGeneUMI_All` resolved to the same variant as
`MultiGeneUMI`, which is neither what STAR does nor what the option is
documented to do. Of the three available behaviours it was the only one nobody
had asked for.
In STAR the option is a no-op: it is parsed and stored, but its consumption site
tests only the `MultiGeneUMI` flag, so selecting it leaves the filter entirely
off. Documented, it removes a UMI seen in more than one gene from *all* of them,
rather than from the losers only.
`UmiFiltering::MultiGeneUmiAll` now exists and does the documented thing: a UMI
appearing in several genes is evidence of a collision or of chimeric
amplification, so it is discarded outright rather than attributed to whichever
gene happened to read deepest. Single-gene UMIs are untouched, which the test
checks across every mode.
Raised upstream as #144 before changing it, since "be faithful to STAR" and "do
what the flag says" genuinely point in opposite directions here.
Also adds `docs-old/dev/divergences.md`, recording this and the homopolymer-UMI
rule, so deliberate differences are written down rather than rediscovered as
surprises in a differential run.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
* docs(changelog): record the MultiGeneUMI_All fix
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
* docs(divergence): correct the MultiGeneUMI_All entry, defer the
homopolymer one
---------
Co-authored-by: Benjamin Demaille <benjamin.demaille@icloud.com>
Co-authored-by: Claude Opus 5 (1M context) <noreply@anthropic.com>
|
Closed by #212 |
Stacked on #174, which is stacked on #173. Merge in that order; this rebases to one commit afterwards.
Finishes the
MultiGeneUMI_CRrule. #173 fixed the tie; this adds the condition I explicitly deferred there, and it needed the restructure I said it would.What STAR actually does
Two things, in this order (
SoloFeature_collapseUMIall.cpp:134-148, then:203-235):UMIs are corrected within each gene first, and the per-gene read totals are recorded twice — once under the raw UMI (
umiGeneMapCount0), once under the corrected one (umiGeneMapCount).Then ownership is decided on the corrected map, subject to two conditions:
The second condition exists because correction moves reads between UMIs. A gene can end up winning only because correction folded a neighbouring UMI onto it, and STAR refuses that win instead of counting it.
Why this could not be a smaller change
The generic path here filters multi-gene UMIs first and corrects afterwards. That order cannot express either condition: by the time correction runs, ownership is already decided. So
MultiGeneUMI_CRnow takes its own path through the per-cell loop, which mirrors STAR — the flag is only valid with--soloUMIdedup 1MM_CR, so no other combination is affected.cellranger_1mm_mapexposes the mappingcellranger_1mmwas already computing and discarding.Measured against CellRanger 10.0.0
Real
cellranger count, not a proxy:cellranger mkrefon the yeast reference, thencellranger counton the fixture from #172. With #165 and #173 also applied:Entries CellRanger has and we do not: 29 → 7.
#165 is a precondition for these numbers and should merge first; without it the same stack sits at 15 439, +2.17%. With it, STAR 2.7.11b on the same fixture with the same flags counts 15 124, so it sits +13 from CellRanger while this sits +5. All three agree to within a fraction of a percent. The full table is in the comment thread below.
Verification
multi_gene_umi_cr_rejects_a_winner_that_only_wins_after_correction— two UMIs one substitution apart, arranged so a gene wins after correction and loses before itmulti_gene_umi_cr_drops_a_tie_entirelyandmulti_gene_umi_cr_keeps_top_genefrom fix(solo): MultiGeneUMI_CR was inert — a tied UMI goes to nobody, not everybody #173 still pass unchangedGate: 562 lib + 26 integration tests,
cargo clippy --all-targets -- -D warnings,cargo fmt --check, all green.What is still different
7 entries CellRanger has that we do not, 19 the other way, of 13 709. I have not chased them and will not guess at a cause here.
The honest limit of this method: anything left is small enough that telling a genuine algorithmic difference from a fixture artefact needs a second, differently-shaped dataset. That belongs with the fixture work in #172 rather than with this rule.